// Solve the Momentum equation

MRF.correctBoundaryVelocity(U);

if (hasFvSource_)
{
    daFvSourcePtr_->calcFvSource(fvSource);
}

tmp<fvVectorMatrix> tUEqn(
    fvm::div(phi, U)
    + MRF.DDt(rho, U)
    + turbulencePtr_->divDevRhoReff(U)
    - fvSource);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();

// get the solver performance info such as initial
// and final residuals
SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));

this->primalResidualControl<vector>(solverU, printToScreen, printInterval, "U");

// bound U
DAUtility::boundVar(allOptions, U, printToScreen);
